Cervezas Artesanas
Debido a mi afición por la Cerveza Artesana, he seleccionado de Kaggle el dataset “beer-recipes” https://www.kaggle.com/jtrofe/beer-recipes
Este es un conjunto de datos de 75.000 cervezas caseras con más de 176 estilos diferentes. Los registros de cerveza son reportados por los usuarios y se clasifican de acuerdo a uno de los 176 estilos diferentes. Estas recetas entran tanto o tan poco en detalle como el usuario proporcionó, pero hay al menos 5 columnas útiles donde se introdujeron datos para cada una: Gravedad Original, Gravedad Final, ABV, IBU, y Color
Datos básicos de las recetas de cerveza presentadas por los usuarios y raspadas de Brewer’s Friend. Todas las columnas están estandarizadas excepto “Método de cebado” y “Cantidad de cebado”, que parece que permite a los usuarios escribir lo que quieran. Agradecimientos
El sitio Brewer’s Friend permite a los usuarios compartir sus recetas de cerveza casera. Este conjunto de datos contiene una selección de las recetas subidas hasta ahora. Inspiración
Sería interesante ver si los datos proporcionados son suficientes para definir cada clase o si hay patrones no descubiertos. En el futuro podría ser posible revisar y raspar información más detallada para cada receta, como la levadura y los lúpulos específicos utilizados.
Alcohol por Volumen. Expresado como porcentaje (%), se utiliza como la medida de la intensidad del alcohol, basado en la proporción de su contenido por volumen de cerveza.
Los estilos de cerveza son las categorías mediante las cuales se identifican y clasifican las cervezas en base a sus características de apariencia, aroma, sabor y sensación en boca, establecidas mediante un rango aproximado de estadísticas vitales.
Unidad Internacional de amargor. Es la medida utilizada usada para definir el amargor de una cerveza. Un IBU equivale a un miligramo de iso-alfa-ácidos disueltos en un litro de cerveza. Cuanto mayor es el valor, más amarga es la cerveza.
El color de la cerveza, depende principalmente, del tipo o tipos de maltas que se utilizan durante su elaboración. Hay 2 metodos el Americano (SRM) y el europeo (EBC)
Analizar si hay relación de los ingredientes y caracteristicas de cada cerveza que podamos obtener un patron de clasificación de las mismas, a pesar de su estilo.
# Limpiamos el entorno de Trabajo
rm(list=ls())
# Limpiamos la consola
cat("\014")# Comprobamos que está bien establecido el directorio
getwd()## [1] "/home/oscar/Documentos/R/Cerveza"dir()## [1] "cerveza-cluster.Rmd" "data" "header.html"#indicamos el directorio de trabajo
setwd("~/Documentos/R/Cerveza")packages <- c("psych", "dplyr", "skimr", "magrittr", "stringr", "epiDisplay", "summarytools","GGally", "ggplot2", "corrplot" , "githubinstall", "e1071", "purrr")
newpack = packages[!(packages %in% installed.packages()[,"Package"])]
if(length(newpack)) install.packages(newpack)
a=lapply(packages, library, character.only=TRUE)Utilizaré un dataset proveniente de “kaggle”
Debido a la restricción que ofrece para poder descargar vía R su Dataset, voy a utilizar la terminal y mediante bash lo voy a descargar y descomprimir.
Este semestre estoy tambien matriculado en “Programación en Scripting”
Para obtener la url y descargar el dataset utilizo la extensión para Chrome wget/curl y realizo la siguientes etapas:
1 - Instalo la extensión.
2 - Empiezo a descargar y lo paro.
3 - Consigo el link de la extensión.
4 - wget o curl + el link.Descargo el fichero: Instalar kaggle - utilizado el terminal mediante comandos bash.
# pip install kaggleObtengo el dataset - utilizado el terminal mediante comandos bash.
Eliminio los posibles ficheros que pudiera tener y descomprimo. - utilizado el terminal mediante comandos bash.
Descomprimo.
unzip(zipfile = "/home/oscar/Descargas/beer-recipes.zip", exdir = "/home/oscar/Documentos/R/Cerveza/data")El dataset que voy a utilizar es el “recipeData.csv” de 13mb. No utilizo el styleData ya que dichos datos están incorporados en el anterior.
recipeData <- read.csv("~/Documentos/R/Cerveza/data/recipeData.csv", encoding = "latin1")
head(recipeData)## BeerID Name
## 1 1 Vanilla Cream Ale
## 2 2 Southern Tier Pumking clone
## 3 3 Zombie Dust Clone - EXTRACT
## 4 4 Zombie Dust Clone - ALL GRAIN
## 5 5 Bakke Brygg Belgisk Blonde 50 l
## 6 6 Sierra Nevada Pale Ale Clone
## URL
## 1 /homebrew/recipe/view/1633/vanilla-cream-ale
## 2 /homebrew/recipe/view/16367/southern-tier-pumking-clone
## 3 /homebrew/recipe/view/5920/zombie-dust-clone-extract
## 4 /homebrew/recipe/view/5916/zombie-dust-clone-all-grain
## 5 /homebrew/recipe/view/89534/bakke-brygg-belgisk-blonde-50-l
## 6 /homebrew/recipe/view/28546/sierra-nevada-pale-ale-clone
## Style StyleID Size.L. OG FG ABV IBU
## 1 Cream Ale 45 21.77 1.055 1.013 5.48 17.65
## 2 Holiday/Winter Special Spiced Beer 85 20.82 1.083 1.021 8.16 60.65
## 3 American IPA 7 18.93 1.063 1.018 5.91 59.25
## 4 American IPA 7 22.71 1.061 1.017 5.80 54.48
## 5 Belgian Blond Ale 20 50.00 1.060 1.010 6.48 17.84
## 6 American Pale Ale 10 24.61 1.055 1.013 5.58 40.12
## Color BoilSize BoilTime BoilGravity Efficiency MashThickness SugarScale
## 1 4.83 28.39 75 1.038 70 N/A Specific Gravity
## 2 15.64 24.61 60 1.07 70 N/A Specific Gravity
## 3 8.98 22.71 60 N/A 70 N/A Specific Gravity
## 4 8.50 26.50 60 N/A 70 N/A Specific Gravity
## 5 4.57 60.00 90 1.05 72 N/A Specific Gravity
## 6 8.00 29.34 70 1.047 79 N/A Specific Gravity
## BrewMethod PitchRate PrimaryTemp PrimingMethod PrimingAmount UserId
## 1 All Grain N/A 17.78 corn sugar 4.5 oz 116
## 2 All Grain N/A N/A N/A N/A 955
## 3 extract N/A N/A N/A N/A NA
## 4 All Grain N/A N/A N/A N/A NA
## 5 All Grain N/A 19 Sukkerlake 6-7 g sukker/l 18325
## 6 All Grain 1 N/A N/A N/A 5889print("Nombre de las columnas: ")## [1] "Nombre de las columnas: "names(recipeData)## [1] "BeerID" "Name" "URL" "Style"
## [5] "StyleID" "Size.L." "OG" "FG"
## [9] "ABV" "IBU" "Color" "BoilSize"
## [13] "BoilTime" "BoilGravity" "Efficiency" "MashThickness"
## [17] "SugarScale" "BrewMethod" "PitchRate" "PrimaryTemp"
## [21] "PrimingMethod" "PrimingAmount" "UserId"print("Clase: ")## [1] "Clase: "class(recipeData)## [1] "data.frame"print("Dimensión: ")## [1] "Dimensión: "dim(recipeData)## [1] 73814 23print("Número de columnas: ")## [1] "Número de columnas: "length(recipeData)## [1] 23print("Estructura: ")## [1] "Estructura: "str(recipeData)## 'data.frame': 73814 obs. of 23 variables:
## $ BeerID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Name : chr "Vanilla Cream Ale" "Southern Tier Pumking clone" "Zombie Dust Clone - EXTRACT" "Zombie Dust Clone - ALL GRAIN" ...
## $ URL : chr "/homebrew/recipe/view/1633/vanilla-cream-ale" "/homebrew/recipe/view/16367/southern-tier-pumking-clone" "/homebrew/recipe/view/5920/zombie-dust-clone-extract" "/homebrew/recipe/view/5916/zombie-dust-clone-all-grain" ...
## $ Style : chr "Cream Ale" "Holiday/Winter Special Spiced Beer" "American IPA" "American IPA" ...
## $ StyleID : int 45 85 7 7 20 10 86 45 129 86 ...
## $ Size.L. : num 21.8 20.8 18.9 22.7 50 ...
## $ OG : num 1.05 1.08 1.06 1.06 1.06 ...
## $ FG : num 1.01 1.02 1.02 1.02 1.01 ...
## $ ABV : num 5.48 8.16 5.91 5.8 6.48 5.58 7.09 5.36 5.77 8.22 ...
## $ IBU : num 17.6 60.6 59.2 54.5 17.8 ...
## $ Color : num 4.83 15.64 8.98 8.5 4.57 ...
## $ BoilSize : num 28.4 24.6 22.7 26.5 60 ...
## $ BoilTime : int 75 60 60 60 90 70 90 75 75 60 ...
## $ BoilGravity : chr "1.038" "1.07" "N/A" "N/A" ...
## $ Efficiency : num 70 70 70 70 72 79 75 70 73 70 ...
## $ MashThickness: chr "N/A" "N/A" "N/A" "N/A" ...
## $ SugarScale : chr "Specific Gravity" "Specific Gravity" "Specific Gravity" "Specific Gravity" ...
## $ BrewMethod : chr "All Grain" "All Grain" "extract" "All Grain" ...
## $ PitchRate : chr "N/A" "N/A" "N/A" "N/A" ...
## $ PrimaryTemp : chr "17.78" "N/A" "N/A" "N/A" ...
## $ PrimingMethod: chr "corn sugar" "N/A" "N/A" "N/A" ...
## $ PrimingAmount: chr "4.5 oz" "N/A" "N/A" "N/A" ...
## $ UserId : int 116 955 NA NA 18325 5889 1051 116 116 NA ...Primera eliminación de columnas no necesarias
data <- dplyr::select(recipeData,-BeerID, -Name, -URL, -PrimingMethod, -PrimingMethod,-UserId, -StyleID)
summary(data)## Style Size.L. OG FG
## Length:73814 Min. : 1.00 Min. : 1.000 Min. :-0.003
## Class :character 1st Qu.: 18.93 1st Qu.: 1.051 1st Qu.: 1.011
## Mode :character Median : 20.82 Median : 1.058 Median : 1.013
## Mean : 43.94 Mean : 1.406 Mean : 1.076
## 3rd Qu.: 23.66 3rd Qu.: 1.069 3rd Qu.: 1.017
## Max. :9200.00 Max. :34.035 Max. :23.425
## NA's :1 NA's :1 NA's :1
## ABV IBU Color BoilSize
## Min. : 0.000 Min. : 0.00 Min. : 0.00 Min. : 1.00
## 1st Qu.: 5.080 1st Qu.: 23.37 1st Qu.: 5.17 1st Qu.: 20.82
## Median : 5.790 Median : 35.77 Median : 8.44 Median : 27.44
## Mean : 6.137 Mean : 44.28 Mean : 13.41 Mean : 49.73
## 3rd Qu.: 6.830 3rd Qu.: 56.39 3rd Qu.: 16.79 3rd Qu.: 30.00
## Max. :54.720 Max. :3409.30 Max. :186.00 Max. :9700.00
## NA's :1 NA's :1 NA's :1 NA's :1
## BoilTime BoilGravity Efficiency MashThickness
## Min. : 0.00 Length:73814 Min. : 0.00 Length:73814
## 1st Qu.: 60.00 Class :character 1st Qu.: 65.00 Class :character
## Median : 60.00 Mode :character Median : 70.00 Mode :character
## Mean : 65.08 Mean : 66.35
## 3rd Qu.: 60.00 3rd Qu.: 75.00
## Max. :240.00 Max. :100.00
## NA's :1 NA's :1
## SugarScale BrewMethod PitchRate PrimaryTemp
## Length:73814 Length:73814 Length:73814 Length:73814
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## PrimingAmount
## Length:73814
## Class :character
## Mode :character
##
##
##
## glimpse(data)## Rows: 73,814
## Columns: 17
## $ Style <chr> "Cream Ale", "Holiday/Winter Special Spiced Beer", "Ame…
## $ Size.L. <dbl> 21.77, 20.82, 18.93, 22.71, 50.00, 24.61, 22.71, 20.82,…
## $ OG <dbl> 1.055, 1.083, 1.063, 1.061, 1.060, 1.055, 1.072, 1.054,…
## $ FG <dbl> 1.013, 1.021, 1.018, 1.017, 1.010, 1.013, 1.018, 1.014,…
## $ ABV <dbl> 5.48, 8.16, 5.91, 5.80, 6.48, 5.58, 7.09, 5.36, 5.77, 8…
## $ IBU <dbl> 17.65, 60.65, 59.25, 54.48, 17.84, 40.12, 268.71, 19.97…
## $ Color <dbl> 4.83, 15.64, 8.98, 8.50, 4.57, 8.00, 6.33, 5.94, 34.76,…
## $ BoilSize <dbl> 28.39, 24.61, 22.71, 26.50, 60.00, 29.34, 30.28, 28.39,…
## $ BoilTime <int> 75, 60, 60, 60, 90, 70, 90, 75, 75, 60, 90, 90, 60, 60,…
## $ BoilGravity <chr> "1.038", "1.07", "N/A", "N/A", "1.05", "1.047", "N/A", …
## $ Efficiency <dbl> 70, 70, 70, 70, 72, 79, 75, 70, 73, 70, 74, 70, 70, 30,…
## $ MashThickness <chr> "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "1.4",…
## $ SugarScale <chr> "Specific Gravity", "Specific Gravity", "Specific Gravi…
## $ BrewMethod <chr> "All Grain", "All Grain", "extract", "All Grain", "All …
## $ PitchRate <chr> "N/A", "N/A", "N/A", "N/A", "N/A", "1", "N/A", "N/A", "…
## $ PrimaryTemp <chr> "17.78", "N/A", "N/A", "N/A", "19", "N/A", "N/A", "N/A"…
## $ PrimingAmount <chr> "4.5 oz", "N/A", "N/A", "N/A", "6-7 g sukker/l", "N/A",…# library(psych)
describe(data)## vars n mean sd median trimmed mad min max
## Style* 1 73814 NaN NA NA NaN NA Inf -Inf
## Size.L. 2 73813 43.94 180.43 20.82 22.33 2.80 1.00 9200.00
## OG 3 73813 1.41 2.20 1.06 1.06 0.01 1.00 34.03
## FG 4 73813 1.08 0.43 1.01 1.01 0.00 0.00 23.42
## ABV 5 73813 6.14 1.88 5.79 5.94 1.22 0.00 54.72
## IBU 6 73813 44.28 42.95 35.77 39.31 21.93 0.00 3409.30
## Color 7 73813 13.41 11.95 8.44 11.11 6.17 0.00 186.00
## BoilSize 8 73813 49.73 193.31 27.44 26.63 6.76 1.00 9700.00
## BoilTime 9 73813 65.08 15.02 60.00 63.37 0.00 0.00 240.00
## BoilGravity* 10 73814 1.35 1.93 1.05 1.05 0.01 0.00 52.60
## Efficiency 11 73813 66.35 14.09 70.00 68.43 7.41 0.00 100.00
## MashThickness* 12 73814 2.13 1.68 1.50 1.95 0.37 0.00 100.00
## SugarScale* 13 73814 NaN NA NA NaN NA Inf -Inf
## BrewMethod* 14 73814 NaN NA NA NaN NA Inf -Inf
## PitchRate* 15 73814 0.75 0.39 0.75 0.70 0.37 0.00 2.00
## PrimaryTemp* 16 73814 19.18 4.22 20.00 19.39 1.65 -17.78 114.00
## PrimingAmount* 17 73813 64.33 133.36 8.00 46.65 10.75 0.00 2500.00
## range skew kurtosis se
## Style* -Inf NA NA NA
## Size.L. 9199.00 15.87 397.42 0.66
## OG 33.03 6.70 47.28 0.01
## FG 23.43 9.78 174.25 0.00
## ABV 54.72 4.95 86.81 0.01
## IBU 3409.30 19.20 1113.91 0.16
## Color 186.00 1.59 2.20 0.04
## BoilSize 9699.00 15.85 402.07 0.71
## BoilTime 240.00 1.38 12.60 0.06
## BoilGravity* 52.60 7.38 67.09 0.01
## Efficiency 100.00 -1.48 1.79 0.05
## MashThickness* 100.00 16.85 540.01 0.01
## SugarScale* -Inf NA NA NA
## BrewMethod* -Inf NA NA NA
## PitchRate* 2.00 1.05 0.74 0.00
## PrimaryTemp* 131.78 3.26 63.79 0.02
## PrimingAmount* 2500.00 11.41 190.89 0.49# valores vacios
print("Mostrar variables con campos na")## [1] "Mostrar variables con campos na"which (is.na(data))## [1] 147628 221442 295256 369070 442884 516698 590512 664326 811954
## [10] 1182422print("Es cierto que hay valores na?")## [1] "Es cierto que hay valores na?"any(is.na(data))## [1] TRUEprint("Suma valores na")## [1] "Suma valores na"sum(is.na(data))## [1] 10print("Mostrar variables con campos na")## [1] "Mostrar variables con campos na"colSums(is.na(data))## Style Size.L. OG FG ABV
## 0 1 1 1 1
## IBU Color BoilSize BoilTime BoilGravity
## 1 1 1 1 0
## Efficiency MashThickness SugarScale BrewMethod PitchRate
## 1 0 0 0 0
## PrimaryTemp PrimingAmount
## 0 1print("Mostrar variables con datos N/A")## [1] "Mostrar variables con datos N/A"colSums(data=="N/A")## Style Size.L. OG FG ABV
## 595 NA NA NA NA
## IBU Color BoilSize BoilTime BoilGravity
## NA NA NA NA 2990
## Efficiency MashThickness SugarScale BrewMethod PitchRate
## NA 29843 0 0 39228
## PrimaryTemp PrimingAmount
## 22646 NAprint("Valores con integrogación")## [1] "Valores con integrogación"colSums(data=="?")## Style Size.L. OG FG ABV
## 0 NA NA NA NA
## IBU Color BoilSize BoilTime BoilGravity
## NA NA NA NA 0
## Efficiency MashThickness SugarScale BrewMethod PitchRate
## NA 0 0 0 0
## PrimaryTemp PrimingAmount
## 0 NAprint("Valores con espacios en blanco")## [1] "Valores con espacios en blanco"colSums(data==" ")## Style Size.L. OG FG ABV
## 0 NA NA NA NA
## IBU Color BoilSize BoilTime BoilGravity
## NA NA NA NA 0
## Efficiency MashThickness SugarScale BrewMethod PitchRate
## NA 0 0 0 0
## PrimaryTemp PrimingAmount
## 0 NAdata$Style <- gsub("N/A", NA, data$Style, fixed = TRUE)
data$BoilSize <- gsub("N/A", NA, data$BoilSize, fixed = TRUE)
data$PitchRate <- gsub("N/A",NA, data$PitchRate, fixed = TRUE)
data$Style <- gsub("N/A",NA, data$Style, fixed = TRUE)
data$BoilGravity <- gsub("N/A",NA, data$BoilGravity, fixed = TRUE)# me quedo con los datos completos
completos <- complete.cases(df)
datos <- df[completos,]# valores vacios
print("Mostrar variables con campos na")## [1] "Mostrar variables con campos na"which (is.na(datos))## integer(0)print("Es cierto que hay valores na?")## [1] "Es cierto que hay valores na?"any(is.na(datos))## [1] FALSEprint("Suma valores na")## [1] "Suma valores na"sum(is.na(datos))## [1] 0print("Mostrar variables con campos na")## [1] "Mostrar variables con campos na"colSums(is.na(datos))## Style Size.L. OG FG ABV IBU
## 0 0 0 0 0 0
## Color BoilSize BoilTime BoilGravity Efficiency SugarScale
## 0 0 0 0 0 0
## BrewMethod
## 0print("Mostrar variables con datos N/A")## [1] "Mostrar variables con datos N/A"colSums(datos=="N/A")## Style Size.L. OG FG ABV IBU
## 0 0 0 0 0 0
## Color BoilSize BoilTime BoilGravity Efficiency SugarScale
## 0 0 0 0 0 0
## BrewMethod
## 0print("Valores con integrogación")## [1] "Valores con integrogación"colSums(datos=="?")## Style Size.L. OG FG ABV IBU
## 0 0 0 0 0 0
## Color BoilSize BoilTime BoilGravity Efficiency SugarScale
## 0 0 0 0 0 0
## BrewMethod
## 0print("Valores con espacios en blanco")## [1] "Valores con espacios en blanco"colSums(datos==" ")## Style Size.L. OG FG ABV IBU
## 0 0 0 0 0 0
## Color BoilSize BoilTime BoilGravity Efficiency SugarScale
## 0 0 0 0 0 0
## BrewMethod
## 0# Modificar los estilos en mayuscula
# library(stringr)
datos$Style <- str_to_upper(datos$Style, locale = "en") #vocablos en ingles
# el resto de vairables character en minusculas
datos$BoilSize <- str_to_lower(datos$BoilSize, locale = "en")
datos$BoilGravity <- str_to_lower(datos$BoilGravity, locale = "en")
# Elimino espacios en blanco
datos$Style <- str_trim(datos$Style, side = "both")
# Separo la columna Style en 2 partes
datos <- datos %>%
tidyr::separate(Style, c("Estilo", "Clase", "Otros"), sep = " ", fill="right")
# Elimino la columna "Otros"
datos <- dplyr::select(datos, -Otros)
# Sustituyo valores na de la columna clase con los valores de la columna Estilo
datos$Clase[is.na(datos$Clase)]<- as.character(datos$Estilo[is.na(datos$Clase)])
# Ref. https://stackoverflow.com/questions/15629885/replace-na-in-column-with-value-in-adjacent-column
# Relleno el espacio en blanco en Clase indicando "otros"
datos$Clase <- sub("^$", "OTHERS", datos$Clase)
# Ref. https://stackoverflow.com/questions/21243588/replace-blank-cells-with-character
# Elimino caracteres especiales
datos$Clase <- gsub("[[:punct:]]", "", datos$Clase)Reducción dataset en 5 tipos de cerverza
# Seleccionamos las Clases de cerveza por nº de frecuencia, hasta que llega a 50% del total y poder obtener clusteres bien definicos
top_data1 <- filter(datos, Clase %in% c("IPA","PALE","STOUT","ALE","PORTER"))Reduccion de observaciones
# Dado que el dataset es de tamaño muy grande, reduzco a 1000 primeras filas aleatorias sin reposición
top_data <- top_data1[sample(nrow(top_data1), 1000, replace=FALSE),]Nuevo resumen de datos
skimr::skim(top_data)| Name | top_data |
| Number of rows | 1000 |
| Number of columns | 14 |
| _______________________ | |
| Column type frequency: | |
| character | 6 |
| numeric | 8 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Estilo | 0 | 1 | 3 | 13 | 0 | 20 | 0 |
| Clase | 0 | 1 | 3 | 6 | 0 | 5 | 0 |
| BoilSize | 0 | 1 | 1 | 7 | 0 | 207 | 0 |
| BoilGravity | 0 | 1 | 2 | 5 | 0 | 135 | 0 |
| SugarScale | 0 | 1 | 5 | 16 | 0 | 2 | 0 |
| BrewMethod | 0 | 1 | 4 | 12 | 0 | 4 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Size.L. | 0 | 1 | 36.81 | 116.43 | 3.03 | 18.93 | 20.82 | 23.00 | 1994.91 | ▇▁▁▁▁ |
| OG | 0 | 1 | 1.30 | 1.76 | 1.02 | 1.05 | 1.06 | 1.07 | 16.54 | ▇▁▁▁▁ |
| FG | 0 | 1 | 1.06 | 0.36 | 1.00 | 1.01 | 1.01 | 1.02 | 6.48 | ▇▁▁▁▁ |
| ABV | 0 | 1 | 6.23 | 1.42 | 2.25 | 5.26 | 6.05 | 6.93 | 15.92 | ▂▇▁▁▁ |
| IBU | 0 | 1 | 60.29 | 62.38 | 0.00 | 33.94 | 49.63 | 72.28 | 1605.83 | ▇▁▁▁▁ |
| Color | 0 | 1 | 14.07 | 12.98 | 2.25 | 5.87 | 7.92 | 14.65 | 50.00 | ▇▁▁▁▁ |
| BoilTime | 0 | 1 | 64.51 | 13.09 | 20.00 | 60.00 | 60.00 | 60.00 | 240.00 | ▇▂▁▁▁ |
| Efficiency | 0 | 1 | 66.45 | 13.80 | 25.00 | 65.00 | 70.00 | 75.00 | 95.00 | ▂▁▂▇▁ |
Convierto atrubutos character a factor y a numeric
# convierto la columna character a factor
top_data$Estilo <- as.factor(top_data$Estilo)
top_data$Clase <- as.factor(top_data$Clase)
top_data$BoilSize <- as.factor(top_data$BoilSize)
top_data$BoilGravity <- as.factor(top_data$BoilGravity)
# convierto variables con mucho números únicos en númericos
top_data$BoilSize <- as.numeric(top_data$BoilSize)
top_data$BoilGravity <- as.numeric(top_data$BoilGravity)Busco valores únicos
unique(top_data$Estilo)## [1] AMERICAN DRY BALTIC IMPERIAL SWEET
## [6] BELGIAN BLONDE CREAM SPECIALTY ROBUST
## [11] CZECH ENGLISH BROWN KELLERBIER: OATMEAL
## [16] DOUBLE OLD IRISH INTERNATIONAL TROPICAL
## 20 Levels: AMERICAN BALTIC BELGIAN BLONDE BROWN CREAM CZECH DOUBLE ... TROPICALunique(top_data$Clase)## [1] IPA PALE STOUT PORTER ALE
## Levels: ALE IPA PALE PORTER STOUTlength(unique(unlist(top_data[c("Estilo")])))## [1] 20length(unique(unlist(top_data[c("Clase")])))## [1] 5agg <- aggregate(Estilo ~ Clase, top_data, function(x) length(unique(x)))
agg## Clase Estilo
## 1 ALE 3
## 2 IPA 5
## 3 PALE 5
## 4 PORTER 5
## 5 STOUT 7Matriz
pairs.panels(top_data, pch=21,main="Gráfico 01.6: Matriz de Dispersión, Histograma y Correlación")
summarytools::freq(top_data$Clase, order = "freq")## Frequencies
## top_data$Clase
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ------------ ------ --------- -------------- --------- --------------
## IPA 495 49.50 49.50 49.50 49.50
## PALE 241 24.10 73.60 24.10 73.60
## STOUT 122 12.20 85.80 12.20 85.80
## ALE 76 7.60 93.40 7.60 93.40
## PORTER 66 6.60 100.00 6.60 100.00
## <NA> 0 0.00 100.00
## Total 1000 100.00 100.00 100.00 100.00# library(epiDisplay)
epiDisplay::summ(top_data)##
## No. of observations = 1000
##
## Var. name obs. mean median s.d. min. max.
## 1 Estilo 1000 5.189 1 6.018 1 20
## 2 Clase 1000 2.663 2 1.114 1 5
## 3 Size.L. 1000 36.81 20.82 116.43 3.03 1994.91
## 4 OG 1000 1.3 1.06 1.76 1.02 16.54
## 5 FG 1000 1.06 1.01 0.36 1 6.48
## 6 ABV 1000 6.23 6.05 1.42 2.25 15.92
## 7 IBU 1000 60.29 49.63 62.38 0 1605.83
## 8 Color 1000 14.07 7.92 12.98 2.25 50
## 9 BoilSize 1000 91.19 95 47.34 1 207
## 10 BoilTime 1000 64.51 60 13.09 20 240
## 11 BoilGravity 1000 38.88 32.5 22.97 1 135
## 12 Efficiency 1000 66.45 70 13.8 25 95
## 13 SugarScale
## 14 BrewMethodtab1(top_data$Clase, sort.group = "decreasing", cum.percent = TRUE)
## top_data$Clase :
## Frequency Percent Cum. percent
## IPA 495 49.5 49.5
## PALE 241 24.1 73.6
## STOUT 122 12.2 85.8
## ALE 76 7.6 93.4
## PORTER 66 6.6 100.0
## Total 1000 100.0 100.0tab1(top_data$Estilo, sort.group = "decreasing", cum.percent = TRUE)
## top_data$Estilo :
## Frequency Percent Cum. percent
## AMERICAN 579 57.9 57.9
## SPECIALTY 65 6.5 64.4
## IMPERIAL 65 6.5 70.9
## BLONDE 48 4.8 75.7
## ENGLISH 34 3.4 79.1
## OATMEAL 27 2.7 81.8
## DOUBLE 27 2.7 84.5
## SWEET 25 2.5 87.0
## BELGIAN 23 2.3 89.3
## CREAM 21 2.1 91.4
## DRY 18 1.8 93.2
## ROBUST 15 1.5 94.7
## INTERNATIONAL 13 1.3 96.0
## BROWN 10 1.0 97.0
## BALTIC 9 0.9 97.9
## CZECH 8 0.8 98.7
## OLD 7 0.7 99.4
## IRISH 3 0.3 99.7
## TROPICAL 2 0.2 99.9
## KELLERBIER: 1 0.1 100.0
## Total 1000 100.0 100.0tab1(top_data$SugarScale, sort.group = "drecreasing", cum.percent = TRUE)
## top_data$SugarScale :
## Frequency Percent Cum. percent
## Plato 18 1.8 1.8
## Specific Gravity 982 98.2 100.0
## Total 1000 100.0 100.0tab1(top_data$BrewMethod, sort.group = "drecreasing", cum.percent = TRUE)
## top_data$BrewMethod :
## Frequency Percent Cum. percent
## All Grain 656 65.6 65.6
## BIAB 183 18.3 83.9
## extract 120 12.0 95.9
## Partial Mash 41 4.1 100.0
## Total 1000 100.0 100.0Guardo el Datamatrix
write.csv(top_data, file = "data/top_data.csv", row.names = FALSE)Convierto el dataframe en Datamatrix
# Utilizando la libreria skimr obtenemos información adicional
dfmatrix <- top_data
cols <- c("Estilo","Clase","BoilSize","BoilGravity","SugarScale","BrewMethod")
# Convertimos lo valores character a factor
for (i in cols){
dfmatrix[,i] <- as.numeric(dfmatrix[,i])
}
skim(dfmatrix)| Name | dfmatrix |
| Number of rows | 1000 |
| Number of columns | 14 |
| _______________________ | |
| Column type frequency: | |
| numeric | 14 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Estilo | 0 | 1 | 5.19 | 6.02 | 1.00 | 1.00 | 1.00 | 10.00 | 20.00 | ▇▁▂▁▁ |
| Clase | 0 | 1 | 2.66 | 1.11 | 1.00 | 2.00 | 2.00 | 3.00 | 5.00 | ▁▇▃▁▂ |
| Size.L. | 0 | 1 | 36.81 | 116.43 | 3.03 | 18.93 | 20.82 | 23.00 | 1994.91 | ▇▁▁▁▁ |
| OG | 0 | 1 | 1.30 | 1.76 | 1.02 | 1.05 | 1.06 | 1.07 | 16.54 | ▇▁▁▁▁ |
| FG | 0 | 1 | 1.06 | 0.36 | 1.00 | 1.01 | 1.01 | 1.02 | 6.48 | ▇▁▁▁▁ |
| ABV | 0 | 1 | 6.23 | 1.42 | 2.25 | 5.26 | 6.05 | 6.93 | 15.92 | ▂▇▁▁▁ |
| IBU | 0 | 1 | 60.29 | 62.38 | 0.00 | 33.94 | 49.63 | 72.28 | 1605.83 | ▇▁▁▁▁ |
| Color | 0 | 1 | 14.07 | 12.98 | 2.25 | 5.87 | 7.92 | 14.65 | 50.00 | ▇▁▁▁▁ |
| BoilSize | 0 | 1 | 91.19 | 47.34 | 1.00 | 68.00 | 95.00 | 108.00 | 207.00 | ▃▆▇▂▂ |
| BoilTime | 0 | 1 | 64.51 | 13.09 | 20.00 | 60.00 | 60.00 | 60.00 | 240.00 | ▇▂▁▁▁ |
| BoilGravity | 0 | 1 | 38.88 | 22.97 | 1.00 | 26.00 | 32.50 | 43.00 | 135.00 | ▅▇▁▁▁ |
| Efficiency | 0 | 1 | 66.45 | 13.80 | 25.00 | 65.00 | 70.00 | 75.00 | 95.00 | ▂▁▂▇▁ |
| SugarScale | 1000 | 0 | NaN | NA | NA | NA | NA | NA | NA | |
| BrewMethod | 1000 | 0 | NaN | NA | NA | NA | NA | NA | NA |
ggpairs(dfmatrix)
# Relación entre variables
correl <- cor(dfmatrix, method = "pearson")
print("Matriz de correlacion de Pearson")## [1] "Matriz de correlacion de Pearson"correl## Estilo Clase Size.L. OG FG
## Estilo 1.000000000 0.237557775 -0.01088419 0.012448984 0.02799722
## Clase 0.237557775 1.000000000 0.01192920 0.026370989 0.05412004
## Size.L. -0.010884189 0.011929199 1.00000000 0.127756217 0.09608860
## OG 0.012448984 0.026370989 0.12775622 1.000000000 0.93286139
## FG 0.027997220 0.054120040 0.09608860 0.932861393 1.00000000
## ABV 0.199401601 -0.003972842 -0.03665281 -0.024556935 -0.02974529
## IBU -0.005751447 -0.136313374 -0.01806301 -0.023188529 -0.02644654
## Color 0.415466693 0.708621343 -0.01167523 0.009396223 0.03070150
## BoilSize 0.055167831 -0.029795743 0.05264564 0.026152204 0.03293212
## BoilTime 0.060180220 0.046733561 0.01739416 0.038584816 0.07311932
## BoilGravity 0.128738983 0.098392378 0.07629586 0.519083528 0.47878075
## Efficiency -0.020434482 -0.026567755 0.15407467 0.103212289 0.09921376
## SugarScale NA NA NA NA NA
## BrewMethod NA NA NA NA NA
## ABV IBU Color BoilSize BoilTime
## Estilo 0.199401601 -0.005751447 0.415466693 0.055167831 0.06018022
## Clase -0.003972842 -0.136313374 0.708621343 -0.029795743 0.04673356
## Size.L. -0.036652813 -0.018063012 -0.011675228 0.052645638 0.01739416
## OG -0.024556935 -0.023188529 0.009396223 0.026152204 0.03858482
## FG -0.029745286 -0.026446541 0.030701499 0.032932118 0.07311932
## ABV 1.000000000 0.258307211 0.263626415 0.023477048 0.21078223
## IBU 0.258307211 1.000000000 -0.024002402 -0.014568568 0.07074333
## Color 0.263626415 -0.024002402 1.000000000 0.008327619 0.07208234
## BoilSize 0.023477048 -0.014568568 0.008327619 1.000000000 0.09745090
## BoilTime 0.210782230 0.070743333 0.072082343 0.097450898 1.00000000
## BoilGravity 0.369819088 0.083403195 0.194743035 -0.126576850 -0.00856133
## Efficiency 0.096397323 -0.018503861 -0.025605275 0.171964409 0.17102686
## SugarScale NA NA NA NA NA
## BrewMethod NA NA NA NA NA
## BoilGravity Efficiency SugarScale BrewMethod
## Estilo 0.12873898 -0.02043448 NA NA
## Clase 0.09839238 -0.02656776 NA NA
## Size.L. 0.07629586 0.15407467 NA NA
## OG 0.51908353 0.10321229 NA NA
## FG 0.47878075 0.09921376 NA NA
## ABV 0.36981909 0.09639732 NA NA
## IBU 0.08340320 -0.01850386 NA NA
## Color 0.19474304 -0.02560527 NA NA
## BoilSize -0.12657685 0.17196441 NA NA
## BoilTime -0.00856133 0.17102686 NA NA
## BoilGravity 1.00000000 -0.23641183 NA NA
## Efficiency -0.23641183 1.00000000 NA NA
## SugarScale NA NA 1 NA
## BrewMethod NA NA NA 1# Buscamos las variables con mayor correlación, hemos realizado una selección entre 0.40 y 0.85
correl_df <- as.data.frame(as.table(correl))
correl_df_a <- subset(correl_df, abs(Freq) > 0.50 & abs(Freq) <0.85)
correl_df_a[order(correl_df_a$Freq),]## Var1 Var2 Freq
## 53 BoilGravity OG 0.5190835
## 144 OG BoilGravity 0.5190835
## 22 Color Clase 0.7086213
## 100 Clase Color 0.7086213# Correlación cla
cor_clase <- cor(dfmatrix[-2], dfmatrix$Clase)
print(cor_clase)## [,1]
## Estilo 0.237557775
## Size.L. 0.011929199
## OG 0.026370989
## FG 0.054120040
## ABV -0.003972842
## IBU -0.136313374
## Color 0.708621343
## BoilSize -0.029795743
## BoilTime 0.046733561
## BoilGravity 0.098392378
## Efficiency -0.026567755
## SugarScale NA
## BrewMethod NAcor_clase1 <- as.data.frame(cor_clase)
cor_clase1[order(abs(cor_clase1$V1)),]## [1] -0.003972842 0.011929199 0.026370989 -0.026567755 -0.029795743
## [6] 0.046733561 0.054120040 0.098392378 -0.136313374 0.237557775
## [11] 0.708621343 NA NASalvo la relación con Estilo, seleccionolas que me parecen más importantes
Selecciono los mejores valores y reduzco el nº me observaciones
# Elegimos las mejores opciones e incorporo IBU y ABV
top_top <- dplyr::select(top_data, Clase, Color, BoilGravity, FG, Efficiency, IBU, ABV)Eliminamos los dataset que no necesitamos
rm(a, agg, data, datos, df, dfmatrix, recipeData, top_data1)Vamos detectar Outliers
par(mfrow=c(2,3))
boxplot(top_top$Color~top_top$Clase)
boxplot(top_top$BoilGravity~top_top$Clase)
boxplot(top_top$FG~top_top$Clase)
boxplot(top_top$Efficiency~top_top$Clase)
boxplot(top_top$IBU~top_top$Clase)
boxplot(top_top$ABV~top_top$Clase)
boxplot(ABV ~ Clase,
data = top_top,
main = "ABV por Clase")
par(mfrow=c(3,2))
boxplot(top_top$Color~top_top$Clase)
boxplot(top_top$BoilGravity~top_top$Clase)
boxplot(top_top$FG~top_top$Clase)
boxplot(top_top$Efficiency~top_top$Clase)
boxplot(top_top$IBU~top_top$Clase)
boxplot(top_top$ABV~top_top$Clase)
# df es el dataFrame que recibimos (ej. activity)
# colNameData es la columna de los datos (ej. "steps")
# colNameBy es la columna por la que trocearemos (ej. "userId")
outliersReplace <- function(df, colNameData, colNameBy){
# creamos una nueva columna llamada igual que colNameData pero con .R
colNameData.R <- paste(colNameData, "R", sep=".")
df[colNameData.R] <- df[colNameData]
# obtenemos los IDs por los que partir el dataframe
IDs <- unique(df[,c(colNameBy)])
for (id in IDs){
data <- df[df[colNameBy] == id, c(colNameData) ]
Q <- quantile(data)
minimo <- Q[1] # valor minimo
Q1 <- Q[2] # primer cuartil
Me <- Q[3] # mediana
Q3 <- Q[4] # tercer cuartil
maximo <- Q[5] # valor maximo
IQR <- Q3 - Q1
lowLimit <- max(minimo, Q1 - 1.7*IQR)
highLimit <- min(maximo, Q3 + 1.7*IQR)
# todos los valores donde colNameBy es igual a id
# y el valor de colNameData es > Q3 + 1.5 * IQR
# lo reemplazamos por la mediana
df[df[colNameBy] == id & df[colNameData] > highLimit, c(colNameData.R)] <- Me
# lo mismo para el umbral inferior
df[df[colNameBy] == id & df[colNameData] < lowLimit, c(colNameData.R)] <- Me
cat(paste("El", colNameBy, id, "la mediana(", colNameData, ") ==", Me, "\n", sep=" " ))
}
df # devolvemos el valor del dataFrame
}top_top <- outliersReplace(top_top,"Color","Clase")## El Clase IPA la mediana( Color ) == 7.66
## El Clase PALE la mediana( Color ) == 6.47
## El Clase STOUT la mediana( Color ) == 38.09
## El Clase PORTER la mediana( Color ) == 33.185
## El Clase ALE la mediana( Color ) == 4.635top_top <- outliersReplace(top_top,"BoilGravity","Clase")## El Clase IPA la mediana( BoilGravity ) == 35
## El Clase PALE la mediana( BoilGravity ) == 28
## El Clase STOUT la mediana( BoilGravity ) == 38
## El Clase PORTER la mediana( BoilGravity ) == 30.5
## El Clase ALE la mediana( BoilGravity ) == 25.5top_top <- outliersReplace(top_top,"FG","Clase")## El Clase IPA la mediana( FG ) == 1.014
## El Clase PALE la mediana( FG ) == 1.011
## El Clase STOUT la mediana( FG ) == 1.017
## El Clase PORTER la mediana( FG ) == 1.0155
## El Clase ALE la mediana( FG ) == 1.011top_top <- outliersReplace(top_top,"Efficiency","Clase")## El Clase IPA la mediana( Efficiency ) == 70
## El Clase PALE la mediana( Efficiency ) == 70
## El Clase STOUT la mediana( Efficiency ) == 70
## El Clase PORTER la mediana( Efficiency ) == 70
## El Clase ALE la mediana( Efficiency ) == 70top_top <- outliersReplace(top_top,"IBU","Clase")## El Clase IPA la mediana( IBU ) == 68.24
## El Clase PALE la mediana( IBU ) == 39.39
## El Clase STOUT la mediana( IBU ) == 37.9
## El Clase PORTER la mediana( IBU ) == 33.62
## El Clase ALE la mediana( IBU ) == 22.5top_top <- outliersReplace(top_top,"ABV","Clase")## El Clase IPA la mediana( ABV ) == 6.57
## El Clase PALE la mediana( ABV ) == 5.34
## El Clase STOUT la mediana( ABV ) == 5.845
## El Clase PORTER la mediana( ABV ) == 6.02
## El Clase ALE la mediana( ABV ) == 5.135Vemos el resultado
par(mfrow = c(2,1)) # para ponerlos uno encima de otro
boxplot(ABV ~ Clase, data = top_top, main = "Sin reemplazo")
boxplot(ABV.R ~ Clase, data = top_top, main = "Con reemplazo")
par(mfrow=c(3,2))
boxplot(top_top$Color.R~top_top$Clase)
boxplot(top_top$BoilGravity.R~top_top$Clase)
boxplot(top_top$FG.R~top_top$Clase)
boxplot(top_top$Efficiency.R~top_top$Clase)
boxplot(top_top$IBU.R~top_top$Clase)
boxplot(top_top$ABV.R~top_top$Clase)
head(top_top)## Clase Color BoilGravity FG Efficiency IBU ABV Color.R
## 7608 IPA 16.72 26 1.01000 70 40.91 5.34 16.72
## 14197 PALE 5.09 32 1.00900 78 44.45 5.12 5.09
## 2761 STOUT 36.42 129 5.33594 66 44.52 4.88 36.42
## 29583 IPA 5.60 37 1.01500 70 71.78 6.27 5.60
## 1880 STOUT 50.00 25 1.00900 50 41.63 4.88 50.00
## 17415 PALE 9.43 31 1.01100 78 38.89 5.82 9.43
## BoilGravity.R FG.R Efficiency.R IBU.R ABV.R
## 7608 26 1.010 70 40.91 5.34
## 14197 32 1.009 78 44.45 5.12
## 2761 38 1.017 66 44.52 4.88
## 29583 37 1.015 70 71.78 6.27
## 1880 25 1.009 50 41.63 4.88
## 17415 31 1.011 78 38.89 5.82top_top <- dplyr::select(top_top, -Color, -BoilGravity, -FG, -Efficiency, -IBU, -ABV)
names(top_top) <- c("Clase", "Color", "BoilGravity", "FG", "Efficiency", "IBU", "ABV")str(top_top)## 'data.frame': 1000 obs. of 7 variables:
## $ Clase : Factor w/ 5 levels "ALE","IPA","PALE",..: 2 3 5 2 5 3 4 2 2 4 ...
## $ Color : num 16.72 5.09 36.42 5.6 50 ...
## $ BoilGravity: num 26 32 38 37 25 31 28 36 34 40 ...
## $ FG : num 1.01 1.01 1.02 1.01 1.01 ...
## $ Efficiency : num 70 78 66 70 50 78 75 68 80 50 ...
## $ IBU : num 40.9 44.5 44.5 71.8 41.6 ...
## $ ABV : num 5.34 5.12 4.88 6.27 4.88 5.82 6.98 6.2 6.12 6.56 ...skimr::skim(top_top)| Name | top_top |
| Number of rows | 1000 |
| Number of columns | 7 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| numeric | 6 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Clase | 0 | 1 | FALSE | 5 | IPA: 495, PAL: 241, STO: 122, ALE: 76 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Color | 0 | 1 | 12.91 | 12.29 | 2.25 | 5.83 | 7.66 | 11.75 | 50.00 | ▇▁▁▁▁ |
| BoilGravity | 0 | 1 | 33.31 | 13.50 | 1.00 | 26.00 | 31.50 | 38.00 | 119.00 | ▂▇▁▁▁ |
| FG | 0 | 1 | 1.01 | 0.00 | 1.00 | 1.01 | 1.01 | 1.02 | 1.03 | ▂▇▃▁▁ |
| Efficiency | 0 | 1 | 71.04 | 6.10 | 50.00 | 70.00 | 70.00 | 75.00 | 90.00 | ▁▂▇▅▁ |
| IBU | 0 | 1 | 52.98 | 26.62 | 0.00 | 33.86 | 46.78 | 68.24 | 149.47 | ▃▇▅▁▁ |
| ABV | 0 | 1 | 6.18 | 1.22 | 3.66 | 5.28 | 6.04 | 6.84 | 11.51 | ▃▇▃▁▁ |
# ¿Con qué variables tendría sentido un proceso de discretización?
apply(top_top,2, function(x) length(unique(x)))## Clase Color BoilGravity FG Efficiency IBU
## 5 666 87 29 41 869
## ABV
## 414# Separamos el objetivo con el resto de objetivo. El objetivo es la clase de cerveza
top.resto <- top_top[,c(2:7)]normalize <- function(x){
return ((x-min(x))/(max(x)-min(x)))
}
top.resto$Color <- normalize(top.resto$Color)
top.resto$BoilGravity <- normalize(top.resto$BoilGravity)
top.resto$FG <- normalize(top.resto$FG)
top.resto$Efficiency <- normalize(top.resto$Efficiency)
top.resto$IBU <- normalize(top.resto$IBU)
top.resto$ABV <- normalize(top.resto$ABV)ggpairs(top.resto, title ="Correlograma con ggpairs")
# Nice visualization of correlations
ggcorr(top.resto, method = c("everything", "pearson"))
# Scatterplot de Variables
textscatter <- function(data, mapping, ...) {
ggplot(data, mapping, ...) + geom_text()
}
ggpairs(
top_top,
title="Scatterplot de Variables",
columns = c(2:7),
mapping=ggplot2::aes(colour = Clase, label = ABV))
lower = list(continuous = textscatter)Scatterplor FG/BoilGravity
textscatter <- function(data, mapping, ...) {
ggplot(data, mapping, ...) + geom_text()
}
ggpairs(
top_top,
title="Scatterplot of Variables",
columns = c(3:4),
mapping=ggplot2::aes(colour = Clase, label = ABV))
lower = list(continuous = textscatter)Scatterplor ABV/IBU
textscatter <- function(data, mapping, ...) {
ggplot(data, mapping, ...) + geom_text()
}
ggpairs(
top_top,
title="Scatterplot of Variables",
columns = c(6,7),
mapping=ggplot2::aes(colour = Clase, label = ABV))
lower = list(continuous = textscatter)rm(correl, correl_df, correl_df_a)Matriz Dispersión, Histograma y Correlación de datos límipios
pairs.panels(top.resto, pch=21,main="Matriz de Dispersión, Histograma y Correlación")
# Seleccionamos los de mayor corrolación
pairs.panels(top.resto[,c("BoilGravity","FG")], pch=21,main="Matriz de Dispersión, Histograma y Correlación - Boilgravity/FG")
pairs.panels(top.resto[,c("IBU","ABV")], pch=21,main="Matriz de Dispersión, Histograma y Correlación - IBU/ABV")
Aplicaremos PCA a las variables continuas y usaremos la variable categórica para visualizar los PCs más tarde. Noten que en el siguiente código aplicamos una transformación de registro a las variables continuas como sugiere [1] y establecemos el centro y la escala. igual a VERDADERO en la llamada a prcomp para estandarizar las variables antes de la aplicación de PCA:
stats.pca <- prcomp(top.resto)
names(stats.pca)## [1] "sdev" "rotation" "center" "scale" "x"# # Aplicar PCA - escala. = VERDADERO es altamente El método de impresión devuelve la desviación estándar de cada uno de los cuatro PCs, y su rotación (o cargas), que son los coeficientes de las combinaciones lineales de las variables continuas.
# print method
print(stats.pca)## Standard deviations (1, .., p=6):
## [1] 0.27246224 0.21774781 0.15415478 0.11938974 0.09904738 0.08399184
##
## Rotation (n x k) = (6 x 6):
## PC1 PC2 PC3 PC4 PC5
## Color -0.92255229 0.17346261 -0.03487347 0.3291019 -0.08309952
## BoilGravity -0.19055890 -0.25768270 0.13296814 -0.2834284 -0.12049066
## FG -0.27054103 -0.27185737 0.08001683 -0.4335778 0.79118054
## Efficiency 0.04384920 -0.03061757 0.96748654 0.2372918 0.02165595
## IBU 0.08397116 -0.71338183 -0.18384660 0.6452162 0.17124509
## ABV -0.17440822 -0.56546693 0.06979142 -0.3882774 -0.56816895
## PC6
## Color -0.04876709
## BoilGravity 0.88586498
## FG -0.18039432
## Efficiency -0.06582728
## IBU 0.06787383
## ABV -0.41398408El método de trazado devuelve un trazado de las variaciones (eje y) asociadas a las PC (eje x). La figura que figura a continuación es útil para decidir cuántas PCs se deben retener para un análisis más detallado. En este simple caso con sólo 4 PCs esto no es una tarea difícil y podemos ver que las dos primeras PCs explican la mayor parte de la variabilidad de los datos.
# plot method
plot(stats.pca$x[, 1], stats.pca$x[, 2], col = top_top$Clase, main = "PCA", xlab = "PC1", ylab = "PC2")
plot(stats.pca, type = "l")
El método de resumen describe la importancia de las PC. La primera fila describe de nuevo la desviación estándar asociada a cada PC. La segunda fila muestra la proporción de la varianza en los datos explicados por cada componente, mientras que la tercera fila describe la proporción acumulativa de la varianza explicada. Podemos ver allí que los tres primeros PCs representan más del {97%} de la varianza de los datos.
# summary method
summary(stats.pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 0.2725 0.2177 0.1542 0.11939 0.09905 0.08399
## Proportion of Variance 0.4205 0.2686 0.1346 0.08074 0.05557 0.03996
## Cumulative Proportion 0.4205 0.6891 0.8237 0.90447 0.96004 1.00000Podemos usar la función de predicción si observamos nuevos datos y queremos predecir los valores de sus PCs. Sólo para ilustrar, imaginemos que las dos últimas filas de los datos del dataset acaban de llegar y queremos ver cuáles son los valores de sus PCs:
# Predict PCs
# Preción PC de las últimas 2 lineas
predict(stats.pca,
newdata=tail(top.resto, 2))## PC1 PC2 PC3 PC4 PC5 PC6
## 10644 -0.9064684 -0.4643252 0.01369397 -0.001264468 -0.32274953 0.09897626
## 32073 0.2077856 0.1651444 -0.01923754 -0.037109444 -0.04355764 0.02706025Grafico 2 primeros PC, que representa 84,9%
library(devtools)
install_github("vqv/ggbiplot")##
##
checking for file ‘/tmp/Rtmp9siJHw/remotes4872332e4ff8/vqv-ggbiplot-7325e88/DESCRIPTION’ ...
✓ checking for file ‘/tmp/Rtmp9siJHw/remotes4872332e4ff8/vqv-ggbiplot-7325e88/DESCRIPTION’
##
─ preparing ‘ggbiplot’:
## checking DESCRIPTION meta-information ...
✓ checking DESCRIPTION meta-information
##
─ checking for LF line-endings in source and make files and shell scripts
##
─ checking for empty or unneeded directories
## ─ looking to see if a ‘data/datalist’ file should be added
## ─ building ‘ggbiplot_0.55.tar.gz’
##
## library(ggbiplot)
g <- ggbiplot(stats.pca, obs.scale = 1, var.scale = 1,
groups = top_top$Clase, ellipse = TRUE,
circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal',
legend.position = 'top')
print(g)
### PCA on caret package
Como mencioné antes, es posible aplicar primero una transformación Box-Cox para corregir la asimetría, centrar y escalar cada variable y luego aplicar PCA en una llamada a la función de preproceso del paquete caret.
require(caret)
trans = preProcess(top_top[,2:7],
method=c("BoxCox", "center",
"scale", "pca"))
PC = predict(trans, top_top[,2:7])Por defecto, la función mantiene sólo los PC necesarios para explicar al menos el 95% de la variabilidad de los datos, pero esto se puede cambiar a través del umbral de argumento.
# Retained PCs
head(PC, 3)## PC1 PC2 PC3 PC4 PC5 PC6
## 7608 0.9360163 -0.8551249 0.04774747 0.78069526 0.6948588 0.0215140
## 14197 1.3561292 0.9348870 -1.12542728 -0.08015398 0.5043589 0.4718226
## 2761 -0.3161090 -2.0251992 0.33484419 0.23048471 0.0604213 1.1212918# Loadings
trans$rotation## PC1 PC2 PC3 PC4 PC5
## Color -0.29895193 -0.68156718 -0.06441744 0.61482816 0.2527764
## BoilGravity -0.51358728 0.03165505 -0.15377421 -0.45811160 0.5424035
## FG -0.48740502 -0.29644631 -0.06401614 -0.25680399 -0.7699636
## Efficiency -0.04901186 0.30947242 -0.90612573 0.26568134 -0.1007660
## IBU -0.35061739 0.54256730 0.37514905 0.52362289 -0.1387185
## ABV -0.53287661 0.23755267 0.07940508 -0.03747406 0.1402209
## PC6
## Color 0.005893036
## BoilGravity 0.455534598
## FG 0.108007286
## Efficiency 0.004523298
## IBU 0.385392784
## ABV -0.795135127library("factoextra")
res.pca <- prcomp(top.resto, scale = TRUE)
fviz_pca_biplot(res.pca, col.ind = top_top$Clase,
palette = "jco", geom = "point")
La dimensión (Dim.) 1 y 2 retuvieron cerca del 61,1% (40% + 21,1%) del total de la información contenida en el conjunto de datos. Los individuos con un perfil similar se agrupan Las variables que están positivamente correlacionadas están en el mismo lado de las parcelas. Las variables que están negativamente correlacionadas están en el lado opuesto de los gráficos.
En primer lugar, el típico PCA de las muestras sería transponer los datos de manera que las muestras sean filas de la matriz de datos.
library(rafalib)
clase <- as.character(top_top$Clase)
group <- as.fumeric(clase)
# transponemos
# x <- t(top.resto)
# Realizamos un PCA
pc <- prcomp(top.resto)
names(pc)## [1] "sdev" "rotation" "center" "scale" "x"plot(pc$x[, 1], pc$x[, 2], col = group, main = "PCA", xlab = "PC1", ylab = "PC2")
Este PCA equivale a realizar la SVD en los datos centrados, donde el centrado ocurre en las columnas. Podemos utilizar la función de barrido para realizar operaciones arbitrarias en las filas y columnas de una matriz. El segundo argumento especifica que queremos operar en las columnas (1 se usaría para las filas), y el tercero y cuarto argumentos especifican que queremos restar la media de las columnas.
cx <- sweep(top.resto, 2, colMeans(top.resto), "-")
sv <- svd(cx)
names(sv)## [1] "d" "u" "v"plot(sv$u[, 1], sv$u[, 2], col = group, main = "SVD", xlab = "U1", ylab = "U2")
Así que las columnas de U de la SVD corresponden a los componentes principales x en el PCA. Además, la matriz V de la SVD equivale a la matriz de rotación devuelta por el PCA.
plot(sv$u[, 1], sv$u[, 2], col = group, main = "SVD", xlab = "U1", ylab = "U2")
plot(pc$x[, 1], pc$x[, 2], col = group, main = "PCA", xlab = "PC1", ylab = "PC2")
plot(sv$u[, 1], sv$u[, 2], col = group, main = "SVD", xlab = "U1", ylab = "U2")
Así que las columnas de U de la SVD corresponden a los componentes principales x en el PCA. Además, la matriz V de la SVD equivale a la matriz de rotación devuelta por el PCA.
sv$v[1:5, 1:5]## [,1] [,2] [,3] [,4] [,5]
## [1,] -0.92255229 0.17346261 -0.03487347 0.3291019 -0.08309952
## [2,] -0.19055890 -0.25768270 0.13296814 -0.2834284 -0.12049066
## [3,] -0.27054103 -0.27185737 0.08001683 -0.4335778 0.79118054
## [4,] 0.04384920 -0.03061757 0.96748654 0.2372918 0.02165595
## [5,] 0.08397116 -0.71338183 -0.18384660 0.6452162 0.17124509pc$rotation[1:5, 1:5]## PC1 PC2 PC3 PC4 PC5
## Color -0.92255229 0.17346261 -0.03487347 0.3291019 -0.08309952
## BoilGravity -0.19055890 -0.25768270 0.13296814 -0.2834284 -0.12049066
## FG -0.27054103 -0.27185737 0.08001683 -0.4335778 0.79118054
## Efficiency 0.04384920 -0.03061757 0.96748654 0.2372918 0.02165595
## IBU 0.08397116 -0.71338183 -0.18384660 0.6452162 0.17124509Los elementos diagonales de D de la SVD son proporcionales a las desviaciones estándar devueltas por PCA. La diferencia es que las desviaciones estándar del prcomp son desviaciones estándar de la muestra (el prcomp devuelve estimaciones no sesgadas de la varianza de la muestra, por lo tanto con la corrección n/(n-1)). Los elementos de D se forman tomando la suma de los cuadrados de los componentes principales pero sin dividirlos por el tamaño de la muestra.
print("SDV")## [1] "SDV"head(sv$d^2)## [1] 74.161437 47.366695 23.739932 14.239655 9.800573 7.047574print("PCA")## [1] "PCA"head(pc$sdev^2)## [1] 0.074235673 0.047414109 0.023763696 0.014253909 0.009810383 0.007054629head(sv$d^2/(ncol(top.resto) - 1))## [1] 14.832287 9.473339 4.747986 2.847931 1.960115 1.409515Dividiendo las variaciones por la suma, obtenemos un gráfico de la relación de variación explicada por cada componente principal.
plot(sv$d^2/sum(sv$d^2), xlim = c(0, 7), type = "b", pch = 16, xlab = "Componentes principales",
ylab = "Varianza explicada")
plot(sv$d^2/sum(sv$d^2), type = "b", pch = 16, xlab = "Componentes principales",
ylab = "Varianza explicada")
Obsérvese que, al no centrar los datos antes de ejecutar el svd, se obtiene un trazado ligeramente diferente:
svNoCenter <- svd(top.resto)
plot(pc$x[, 1], pc$x[, 2], col = group, main = "PCA", xlab = "PC1", ylab = "PC2")
points(0, 0, pch = 3, cex = 4, lwd = 4)
plot(svNoCenter$u[, 1], svNoCenter$u[, 2], col = group, main = "SVD no centrada",
xlab = "U1", ylab = "U2")
result<- kmeans(top.resto,4) #aplico k-means con n.. de centroides(k)=4
result$size # Obtenemos el nº de grabaciones por cada cluster## [1] 482 334 60 124result$centers # obtenemos el valor del centro del cluster(3 centers for k=4)## Color BoilGravity FG Efficiency IBU ABV
## 1 0.09684321 0.2247873 0.2866417 0.5266338 0.2630878 0.2380659
## 2 0.12481362 0.3116563 0.3952096 0.5324850 0.5413723 0.4333384
## 3 0.85453054 0.4862288 0.5956989 0.4962500 0.3110390 0.5362208
## 4 0.67371559 0.2594314 0.3694069 0.5211694 0.2272233 0.2354017# Da el vector del cúmulo mostrando el cúmulo donde cae cada registro
head(result$cluster, 100)## 7608 14197 2761 29583 1880 17415 32586 9943 23649 15226 20338 3071 21980
## 1 1 4 2 4 1 4 1 1 4 4 1 1
## 25647 28886 32395 27374 15655 17798 3795 28675 5658 24780 13935 2355 25867
## 4 1 1 1 2 1 1 2 4 1 2 2 2
## 17291 2025 4089 9282 11598 11841 25085 22616 34038 16016 13827 30786 31403
## 1 2 2 3 1 2 4 2 2 1 4 3 1
## 16640 7310 11614 3626 3630 21340 2475 31262 12355 23463 14871 32235 33416
## 2 1 2 4 4 3 1 2 1 2 2 3 1
## 34243 11665 19445 19110 31182 11083 13302 28349 29184 7381 22866 14257 30127
## 1 1 1 4 1 4 1 1 1 2 1 4 1
## 6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338
## 2 2 1 2 1 1 1 4 4 2 1 1 1
## 34258 193 9535 11829 34074 24806 9389 32397 26466 138 32918 4230 32518
## 1 2 2 1 1 1 2 3 2 2 1 3 1
## 30764 9581 31086 27054 53 12057 14422 12010 3190
## 2 1 1 2 2 1 2 2 2Verificamos el resultado del clustering
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,2)], col=result$cluster)
plot(top.resto[c(1,2)], col=top_top$Clase)
plot(top.resto[c(3,4)], col=result$cluster)
plot(top.resto[c(3,4)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(5,6)], col=result$cluster)
plot(top.resto[c(5,6)], col=top_top$Clase)
plot(top.resto[c(1,6)], col=result$cluster)
plot(top.resto[c(1,6)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,3)], col=result$cluster)
plot(top.resto[c(1,3)], col=top_top$Clase)
plot(top.resto[c(2,4)], col=result$cluster)
plot(top.resto[c(2,4)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,5)], col=result$cluster)
plot(top.resto[c(1,5)], col=top_top$Clase)
plot(top.resto[c(3,5)], col=result$cluster)
plot(top.resto[c(3,5)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,4)], col=result$cluster)
plot(top.resto[c(1,4)], col=top_top$Clase)
plot(top.resto[c(3,2)], col=result$cluster)
plot(top.resto[c(3,2)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(6,2)], col=result$cluster)
plot(top.resto[c(6,2)], col=top_top$Clase)
plot(top.resto[c(5,4)], col=result$cluster)
plot(top.resto[c(5,4)], col=top_top$Clase)
table(result$cluster,top_top$Clase)##
## ALE IPA PALE PORTER STOUT
## 1 76 169 234 3 0
## 2 0 326 7 0 1
## 3 0 0 0 10 50
## 4 0 0 0 53 71Según la tabla, el cluster 1 estan IPA, en la 2, PORTER y STOUT en la 3 ninguna, y en la cuarta ALE y PALE
El número total de instancias correctamente clasificadas son: 786 El número total de casos clasificados incorrectamente son: 214 Exactitud = 786/(1000) = 0.786 i.e. nuestro modelo ha alcanzado un 78,6% de exactitud!
result<- kmeans(top.resto,3) #aplico k-means con n.. de centroides(k)=3
result$size # Obtenemos el nº de grabaciones por cada cluster## [1] 60 126 814result$centers # obtenemos el valor del centro del cluster(3 centers for k=3)## Color BoilGravity FG Efficiency IBU ABV
## 1 0.8545305 0.4862288 0.5956989 0.4962500 0.3110390 0.5362208
## 2 0.6687941 0.2569276 0.3684076 0.5228175 0.2255414 0.2336467
## 3 0.1076644 0.2607338 0.3311405 0.5287930 0.3776218 0.3184682# Da el vector del cúmulo mostrando el cúmulo donde cae cada registro
head(result$cluster, 100)## 7608 14197 2761 29583 1880 17415 32586 9943 23649 15226 20338 3071 21980
## 3 3 2 3 2 3 2 3 3 2 2 3 3
## 25647 28886 32395 27374 15655 17798 3795 28675 5658 24780 13935 2355 25867
## 2 3 3 3 3 3 3 3 2 3 3 3 3
## 17291 2025 4089 9282 11598 11841 25085 22616 34038 16016 13827 30786 31403
## 3 3 3 1 3 3 2 3 3 3 2 1 3
## 16640 7310 11614 3626 3630 21340 2475 31262 12355 23463 14871 32235 33416
## 3 3 3 2 2 1 3 3 3 3 3 1 3
## 34243 11665 19445 19110 31182 11083 13302 28349 29184 7381 22866 14257 30127
## 3 3 3 2 3 2 3 3 3 3 3 2 3
## 6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338
## 3 3 3 3 3 3 3 2 2 3 3 3 3
## 34258 193 9535 11829 34074 24806 9389 32397 26466 138 32918 4230 32518
## 3 3 3 3 3 3 3 1 3 3 3 1 3
## 30764 9581 31086 27054 53 12057 14422 12010 3190
## 3 3 3 3 3 3 3 3 3Verificamos el resultado del clustering
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,2)], col=result$cluster)
plot(top.resto[c(1,2)], col=top_top$Clase)
plot(top.resto[c(3,4)], col=result$cluster)
plot(top.resto[c(3,4)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(5,6)], col=result$cluster)
plot(top.resto[c(5,6)], col=top_top$Clase)
plot(top.resto[c(1,6)], col=result$cluster)
plot(top.resto[c(1,6)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,3)], col=result$cluster)
plot(top.resto[c(1,3)], col=top_top$Clase)
plot(top.resto[c(2,4)], col=result$cluster)
plot(top.resto[c(2,4)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,5)], col=result$cluster)
plot(top.resto[c(1,5)], col=top_top$Clase)
plot(top.resto[c(3,5)], col=result$cluster)
plot(top.resto[c(3,5)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,4)], col=result$cluster)
plot(top.resto[c(1,4)], col=top_top$Clase)
plot(top.resto[c(3,2)], col=result$cluster)
plot(top.resto[c(3,2)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(6,2)], col=result$cluster)
plot(top.resto[c(6,2)], col=top_top$Clase)
plot(top.resto[c(5,4)], col=result$cluster)
plot(top.resto[c(5,4)], col=top_top$Clase)
table(result$cluster,top_top$Clase)##
## ALE IPA PALE PORTER STOUT
## 1 0 0 0 10 50
## 2 0 0 0 55 71
## 3 76 495 241 1 1Según la tabla, el cluster 1 estan ALE y PALE, en la 2, claramente IPA y en la tercera, PORTER y STOUT.
El número total de instancias correctamente clasificadas son: 835 El número total de casos clasificados incorrectamente son: 165 Exactitud = 835/(835+165) = 0.835 i.e. nuestro modelo ha alcanzado un 83,5% de exactitud!
result<- kmeans(top.resto,2) #aplico k-means con n.. de centroides(k)=2
result$size # Obtenemos el nº de grabaciones por cada cluster## [1] 183 817result$centers # obtenemos el valor del centro del cluster(3 centers for k=2)## Color BoilGravity FG Efficiency IBU ABV
## 1 0.7343862 0.3340048 0.4438569 0.5131148 0.2545376 0.3342940
## 2 0.1086733 0.2602950 0.3310696 0.5289933 0.3768474 0.3178336# Da el vector del cúmulo mostrando el cúmulo donde cae cada registro
head(result$cluster, 100)## 7608 14197 2761 29583 1880 17415 32586 9943 23649 15226 20338 3071 21980
## 2 2 1 2 1 2 1 2 2 1 1 2 2
## 25647 28886 32395 27374 15655 17798 3795 28675 5658 24780 13935 2355 25867
## 1 2 2 2 2 2 2 2 1 2 2 2 2
## 17291 2025 4089 9282 11598 11841 25085 22616 34038 16016 13827 30786 31403
## 2 2 2 1 2 2 1 2 2 2 1 1 2
## 16640 7310 11614 3626 3630 21340 2475 31262 12355 23463 14871 32235 33416
## 2 2 2 1 1 1 2 2 2 2 2 1 2
## 34243 11665 19445 19110 31182 11083 13302 28349 29184 7381 22866 14257 30127
## 2 2 2 1 2 1 2 2 2 2 2 1 2
## 6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338
## 2 2 2 2 2 2 2 1 1 2 2 2 2
## 34258 193 9535 11829 34074 24806 9389 32397 26466 138 32918 4230 32518
## 2 2 2 2 2 2 2 1 2 2 2 1 2
## 30764 9581 31086 27054 53 12057 14422 12010 3190
## 2 2 2 2 2 2 2 2 2Verificamos el resultado del clustering
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,2)], col=result$cluster)
plot(top.resto[c(1,2)], col=top_top$Clase)
plot(top.resto[c(3,4)], col=result$cluster)
plot(top.resto[c(3,4)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(5,6)], col=result$cluster)
plot(top.resto[c(5,6)], col=top_top$Clase)
plot(top.resto[c(1,6)], col=result$cluster)
plot(top.resto[c(1,6)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,3)], col=result$cluster)
plot(top.resto[c(1,3)], col=top_top$Clase)
plot(top.resto[c(2,4)], col=result$cluster)
plot(top.resto[c(2,4)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,5)], col=result$cluster)
plot(top.resto[c(1,5)], col=top_top$Clase)
plot(top.resto[c(3,5)], col=result$cluster)
plot(top.resto[c(3,5)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,4)], col=result$cluster)
plot(top.resto[c(1,4)], col=top_top$Clase)
plot(top.resto[c(3,2)], col=result$cluster)
plot(top.resto[c(3,2)], col=top_top$Clase)
par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(6,2)], col=result$cluster)
plot(top.resto[c(6,2)], col=top_top$Clase)
plot(top.resto[c(5,4)], col=result$cluster)
plot(top.resto[c(5,4)], col=top_top$Clase)
table(result$cluster,top_top$Clase)##
## ALE IPA PALE PORTER STOUT
## 1 0 0 0 63 120
## 2 76 495 241 3 2Según la tabla, el cluster 2 estan ALE, IPA Y PALE y el cluster 1 están PORTER Y STOUT
El número total de instancias correctamente clasificadas son: 995 El número total de casos clasificados incorrectamente son: 5 Exactitud = 995/(1000) = 0.995 i.e. nuestro modelo ha alcanzado un 99,% de exactitud!
library(factoextra)
fviz_nbclust(x = top.resto, FUNcluster = kmeans, method = "wss", k.max = 7,
diss = get_dist(top.resto, method = "euclidean"), nstart = 30)
# Este mismo análisis también puede realizarse sin recurrir a la función fviz_nbclust().
calcular_totwithinss <- function(n_clusters, datos, iter.max=1000, nstart=50){
# Esta función aplica el algoritmo kmeans y devuelve la suma total de
# cuadrados internos.
cluster_kmeans <- kmeans(centers = n_clusters, x = top.resto, iter.max = iter.max,
nstart = nstart)
return(cluster_kmeans$tot.withinss)
}
# library(purrr)
# Se aplica esta función con para diferentes valores de k
total_withinss <- map_dbl(.x = 1:7,
.f = calcular_totwithinss,
datos = top.resto)
total_withinss## [1] 176.35587 112.79074 85.87822 76.58989 69.04624 64.00273 60.15495data.frame(n_clusters = 1:7, suma_cuadrados_internos = total_withinss) %>%
ggplot(aes(x = n_clusters, y = suma_cuadrados_internos)) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = 1:7) +
labs(title = "Evolución de la suma total de cuadrados intra-cluster") +
theme_bw()
library(cluster)
library(factoextra)
fviz_nbclust(x = top.resto, FUNcluster = pam, method = "wss", k.max = 5,
diss = dist(top.resto, method = "manhattan"))
set.seed(123)
pam_clusters <- pam(x = top.resto, k = 3, metric = "manhattan")
pam_clusters## Medoids:
## ID Color BoilGravity FG Efficiency IBU ABV
## 26633 384 0.04732984 0.2076271 0.2903226 0.5 0.2466716 0.2127389
## 7922 456 0.79057592 0.3135593 0.4193548 0.5 0.2511541 0.3617834
## 8841 727 0.10198953 0.3050847 0.3870968 0.5 0.4388841 0.4433121
## Clustering vector:
## 7608 14197 2761 29583 1880 17415 32586 9943 23649 15226 20338 3071 21980
## 1 1 2 3 2 1 2 3 3 2 2 1 3
## 25647 28886 32395 27374 15655 17798 3795 28675 5658 24780 13935 2355 25867
## 2 1 1 1 3 1 1 3 1 1 3 3 3
## 17291 2025 4089 9282 11598 11841 25085 22616 34038 16016 13827 30786 31403
## 1 3 3 2 1 3 2 3 3 1 2 2 1
## 16640 7310 11614 3626 3630 21340 2475 31262 12355 23463 14871 32235 33416
## 3 1 3 2 2 2 1 3 1 3 3 2 1
## 34243 11665 19445 19110 31182 11083 13302 28349 29184 7381 22866 14257 30127
## 1 1 1 1 1 2 1 1 3 3 3 2 1
## 6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338
## 3 3 1 3 1 1 1 1 2 3 1 3 3
## 34258 193 9535 11829 34074 24806 9389 32397 26466 138 32918 4230 32518
## 1 3 3 1 1 1 1 2 3 3 1 2 1
## 30764 9581 31086 27054 53 12057 14422 12010 3190 32889 28666 33777 5927
## 3 3 1 3 3 1 3 3 3 2 3 1 1
## 6298 21656 3586 28856 31152 19089 13566 33631 29572 12154 26093 25175 26830
## 3 1 3 3 3 3 3 1 3 3 3 3 1
## 11433 4735 23258 734 4014 28804 18054 12340 34510 12765 15247 20063 4216
## 1 1 3 3 2 3 1 2 2 1 1 3 3
## 12138 14321 31450 31817 21195 1331 25300 29627 17785 24632 9379 17188 18878
## 1 3 1 3 2 3 1 3 1 1 1 3 2
## 11621 20658 11331 5421 27666 31845 672 34152 6760 24223 26155 2771 24749
## 3 1 1 3 2 2 3 1 1 1 3 1 1
## 23819 32661 6924 9848 241 19258 2578 28316 5993 29561 10786 27129 9976
## 1 3 1 3 3 3 3 1 3 2 3 3 2
## 18389 12964 2743 20800 12835 4775 4497 31138 16385 10368 17751 13497 33188
## 3 1 2 1 2 1 2 2 3 1 1 3 1
## 449 31363 33331 20086 34142 25046 23742 3112 7969 14044 27219 21428 8305
## 3 3 3 3 3 1 1 3 3 3 1 2 3
## 30962 18861 3222 27706 16803 12674 10451 30355 12300 749 26959 22037 22332
## 1 1 1 1 2 2 2 3 3 1 2 2 3
## 24527 11421 32385 21988 8586 24383 3988 6144 32093 25432 30705 23910 24156
## 3 2 1 2 3 3 3 3 3 3 2 3 1
## 30235 12158 5119 22779 21311 5409 31292 23127 14005 7136 32610 33622 5753
## 1 2 1 1 1 3 2 3 1 1 3 3 1
## 11095 2219 24759 2167 33920 5987 27649 13207 15431 9078 4292 29265 18837
## 3 2 3 3 1 3 3 2 1 1 1 2 1
## 34205 16154 5650 2359 20160 11844 4684 2413 7702 7753 25268 33966 20429
## 2 1 3 2 1 3 1 1 1 3 1 3 1
## 6966 3784 11101 31514 29685 20459 3996 13748 2712 799 5682 17454 33765
## 3 1 1 1 1 1 1 1 3 2 3 3 3
## 2927 30444 6849 24469 10827 23945 23094 30561 31923 16407 32378 8860 27048
## 1 1 1 2 1 1 3 1 1 1 1 1 3
## 12301 31516 32134 14358 21520 16575 28474 18829 6881 18428 29492 23394 22384
## 2 2 1 3 1 1 3 3 1 1 2 1 2
## 16566 34144 9663 22290 23107 17819 9574 31926 4287 17859 3117 15663 21102
## 1 2 2 3 3 3 2 3 3 2 1 1 3
## 8852 34361 2239 9881 32863 13253 23360 9605 28504 6177 28472 2739 26559
## 2 3 1 3 1 3 3 1 3 3 3 3 3
## 34249 29130 14908 7097 32348 7342 7168 5087 24425 912 4856 10280 16888
## 1 3 3 1 1 3 3 2 3 3 1 3 1
## 3673 28685 33271 2857 4669 10363 17876 11955 344 27611 10323 33720 12092
## 1 2 1 2 2 1 3 3 3 1 3 1 2
## 6961 31781 25176 13739 27654 33210 16625 6860 20316 12988 18670 229 32474
## 1 3 1 3 2 3 3 3 3 2 1 1 1
## 15745 17705 31522 27738 9075 8696 8437 32364 12408 30125 25362 25771 10373
## 3 1 1 3 2 3 2 1 3 3 1 3 3
## 18151 19147 25206 1586 33320 6839 26633 27185 1982 31217 2546 25673 6753
## 1 1 1 3 2 1 1 1 1 1 1 1 2
## 1516 33899 21318 1650 22326 16735 34328 18300 31964 18269 34259 34500 23192
## 1 1 2 3 1 3 1 2 1 1 1 1 1
## 18554 18599 8941 31283 20189 10184 16456 16760 19687 28208 17521 9814 14886
## 1 1 3 2 1 1 1 1 3 3 2 1 2
## 12202 3524 18456 10576 30222 32458 8788 3771 13870 28041 15369 21909 26111
## 2 3 3 3 3 3 1 1 3 2 1 3 3
## 13012 24199 20178 18819 20663 24672 18284 10767 16083 6724 16059 2213 31050
## 2 1 1 1 1 3 3 3 1 2 3 3 2
## 20518 27730 3258 12779 3893 17437 8611 808 32463 5017 33636 4679 33911
## 3 3 2 1 1 2 2 2 3 1 3 1 3
## 7922 10060 6767 16630 32680 26584 27149 17364 26270 5337 25240 11554 23935
## 2 1 3 1 1 3 1 3 1 1 1 3 3
## 26947 9753 25245 23613 20092 30595 34651 20098 25328 8088 13914 26013 113
## 1 3 1 3 1 1 2 1 3 3 1 3 1
## 9368 25151 31293 25800 10883 30560 28006 10770 20439 27880 28945 26846 7910
## 1 3 1 1 1 1 3 3 1 3 3 1 1
## 7971 9707 2941 27802 18831 14247 26932 18590 27759 1106 8285 14588 5787
## 3 1 2 3 3 1 3 2 1 1 2 3 2
## 11822 19806 8518 31755 21095 6117 9894 25544 2577 2937 23963 25230 21984
## 2 2 2 3 3 3 2 1 1 3 1 1 1
## 22917 15149 648 26141 12425 11320 24764 25540 22791 27367 1934 10054 20828
## 1 1 1 1 1 1 3 1 1 1 1 1 1
## 24978 20561 4180 30896 34345 32503 1512 11736 16700 9222 33104 23793 5656
## 3 3 3 3 3 2 1 3 3 2 3 3 1
## 19978 12857 2586 31307 22622 9926 33146 19340 9327 20583 34076 17956 3155
## 3 1 3 3 1 3 3 1 3 1 1 3 1
## 11329 19165 27103 8261 27224 26036 24483 19154 30287 22044 31647 11449 5854
## 3 3 2 3 1 3 2 1 1 1 2 1 3
## 29166 13391 26777 16044 13230 3601 12120 19612 6920 13470 16419 29386 21247
## 3 2 1 3 1 3 3 1 3 1 1 3 1
## 22963 2407 20824 13588 30858 20977 10299 31482 7709 11260 15679 23758 11819
## 1 1 1 1 3 3 2 2 1 2 1 1 3
## 25045 25691 34045 14862 6018 25379 1631 9545 17387 6173 23761 22724 6163
## 1 3 1 1 1 1 3 3 3 1 2 1 1
## 12446 20853 24424 23673 19491 23422 11028 8215 28104 8861 18763 30570 17361
## 1 3 1 1 3 3 3 1 1 3 3 3 2
## 32374 23374 20132 32399 15878 445 16729 5732 412 1499 17235 14100 25947
## 1 1 3 3 1 2 1 2 2 3 1 1 2
## 10292 7000 5629 16379 22638 13612 20462 25533 570 29024 3919 30728 32787
## 1 3 1 3 3 3 1 1 3 3 1 2 1
## 23410 27516 25380 12502 25604 18898 3324 8595 19931 7447 23501 24258 2042
## 3 1 3 1 1 1 2 3 1 2 1 3 3
## 3682 22571 11644 26723 22669 9939 27331 21300 23171 15569 27262 28941 16608
## 1 3 3 3 3 3 1 1 1 3 3 1 3
## 9562 16025 32026 21892 19136 30941 3664 3678 11199 11733 2095 32237 28360
## 3 2 3 1 3 1 2 3 1 3 2 1 2
## 32638 17266 2119 14224 21716 6504 7012 34431 17181 5666 5932 31973 20534
## 1 2 3 1 3 3 3 1 2 2 1 1 2
## 23046 32488 21125 931 33409 11610 25964 20894 15437 1391 6637 32265 27052
## 3 3 3 3 1 1 1 1 1 1 2 3 3
## 26027 10040 26016 4510 34448 21468 7739 21162 3162 19627 26508 8841 21940
## 1 1 1 3 3 2 3 1 3 3 3 3 3
## 16547 513 25489 17103 11429 19030 11810 21505 9998 29860 9452 17877 27508
## 3 1 1 3 1 1 1 3 1 1 1 3 1
## 21828 12626 5619 3097 16031 6711 4447 1036 34098 24406 10664 21200 34445
## 1 1 2 3 1 3 1 3 1 3 1 1 1
## 32075 32165 23704 24766 3070 34217 8925 28847 28553 15118 34194 4695 5513
## 2 1 3 1 1 1 3 2 3 3 2 1 1
## 13166 12588 9357 34491 23676 25283 15864 11864 7275 2465 31931 29533 19697
## 3 3 3 1 2 3 2 3 3 1 1 1 1
## 25759 16995 17366 30472 22184 12955 16064 1695 29364 6028 916 7557 32333
## 2 1 1 3 3 3 1 3 3 1 3 3 3
## 31783 31452 21752 20120 23125 33530 1110 25730 15059 19844 1553 32371 17588
## 3 2 2 1 2 1 2 1 3 1 3 1 1
## 4538 32202 14249 17447 25373 5766 911 526 27038 14905 28382 28460 27889
## 3 1 2 2 1 3 3 3 1 1 3 3 2
## 11232 4830 2874 11386 2582 34027 14743 29493 20137 7468 30676 12264 23159
## 1 1 3 3 3 2 1 2 1 1 3 2 1
## 15354 24828 29744 27474 24445 27679 9669 24582 10418 10283 7766 3634 12453
## 3 3 1 1 1 2 1 1 3 1 3 1 2
## 30511 33908 18801 14672 30276 15849 13151 16768 17487 5001 26145 3 19837
## 1 1 1 3 1 3 2 2 1 1 1 1 3
## 15914 14991 4345 30800 27642 24956 20970 12007 7906 11214 30295 15281 17436
## 1 3 3 1 3 2 1 3 3 1 1 1 3
## 31227 27582 10670 32496 5663 17241 3780 11357 34009 19011 5978 8703 823
## 3 3 2 1 3 2 1 1 3 3 3 1 2
## 14090 10016 3186 12099 17328 10395 15366 17275 28546 27587 125 5287 11253
## 3 3 3 1 2 3 1 1 3 1 1 1 1
## 26709 281 31011 32774 32554 34455 25886 3160 1733 24047 6794 15927 7337
## 3 1 3 1 2 3 2 3 3 1 3 2 1
## 14665 33114 14463 8330 30046 31651 34399 6738 17823 9090 15352 2698 552
## 1 3 1 1 1 3 1 1 3 3 1 1 3
## 18355 17726 32563 1155 20706 33276 26471 11263 14454 22573 27896 29512 28705
## 3 3 2 2 3 1 2 1 1 1 1 1 3
## 32753 10429 6939 34636 1671 16005 11891 2984 12058 19912 13816 8441 19128
## 1 2 3 3 3 2 1 3 1 3 3 2 3
## 13181 30698 31294 10625 14264 22777 20290 23830 30499 31949 5488 22444 17920
## 3 1 1 3 1 3 1 1 3 3 1 3 3
## 14478 5895 15246 8264 8056 22496 14190 26273 31471 9079 5768 31649 18404
## 3 2 3 3 3 3 3 1 1 1 1 2 3
## 20520 32982 16987 18360 33832 7999 1628 30171 11393 5456 7231 13510 26937
## 3 1 1 2 3 3 1 1 1 3 1 1 1
## 9155 21569 13945 10662 17172 14920 13940 11413 34070 10387 10644 32073
## 2 3 1 3 1 1 1 3 3 2 2 1
## Objective function:
## build swap
## 0.5337826 0.5104537
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"fviz_cluster(object = pam_clusters, data = top.resto, ellipse.type = "t",
repel = TRUE) +
theme_bw() +
labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
## CLARA
library(cluster)
library(factoextra)
clara_clusters <- clara(x = top.resto, k = 3, metric = "manhattan", stand = TRUE,
samples = 10, pamLike = TRUE)
clara_clusters## Call: clara(x = top.resto, k = 3, metric = "manhattan", stand = TRUE, samples = 10, pamLike = TRUE)
## Medoids:
## Color BoilGravity FG Efficiency IBU ABV
## 32073 0.05089005 0.2288136 0.2258065 0.500 0.2280725 0.2178344
## 12264 0.65005236 0.2500000 0.4032258 0.500 0.2249281 0.3006369
## 31227 0.07790576 0.2966102 0.3548387 0.575 0.4129926 0.3324841
## Objective function: 4.546298
## Clustering vector: Named int [1:1000] 1 1 2 3 2 1 2 3 3 2 2 1 3 2 1 1 3 3 ...
## - attr(*, "names")= chr [1:1000] "7608" "14197" "2761" "29583" "1880" "17415" "32586" ...
## Cluster sizes: 352 176 472
## Best sample:
## [1] 3071 13935 25085 25478 9535 12138 20658 672 12964 17751 12674 6177
## [13] 34249 33271 12408 25673 18269 28208 13012 33911 23613 20092 25151 26846
## [25] 9894 33146 13588 17387 32399 32787 3324 3682 21828 34445 34217 34194
## [37] 12955 31783 12264 30511 31227 5663 12058 13816 11413 32073
##
## Available components:
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"fviz_cluster(object = clara_clusters, ellipse.type = "t", geom = "point",
pointsize = 2.5) +
theme_bw() +
labs(title = "Resultados clustering CLARA") +
theme(legend.position = "none")
## Hierarchical K-means clusteriong
library(factoextra)
# Se obtiene el dendrograma de hierarchical clustering para elegir el número de
# clusters.
set.seed(101)
hc_euclidea_completo <- hclust(d = dist(x = top.resto, method = "euclidean"),
method = "complete")
fviz_dend(x = hc_euclidea_completo, cex = 0.5, main = "Linkage completo",
sub = "Distancia euclídea") +
theme(plot.title = element_text(hjust = 0.5, size = 15))
https://www.r-bloggers.com/experimentation-with-unsupervised-learning/
https://www.r-bloggers.com/computing-and-visualizing-pca-in-r/
https://genomicsclass.github.io/book/pages/pca_svd.html
https://rpubs.com/Nitika/kmeans_Iris
https://rpubs.com/Joaquin_AR/310338
https://www.adictosaltrabajo.com/2019/11/28/deteccion-y-reemplazo-de-outliers-con-r/